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We propose an extension to the Davydov D 2 Ansatz in the dynamics study of the Holstein molec¬ 
ular crystal model with diagonal and off-diagonal exciton-phonon coupling using the Dirac-Frenkel 
time-dependent variational principle. The new trial state by the name of the “multi-D 2 Ansatz is a 
linear combination of Davydov D 2 trial states, and its validity is carefully examined by quantifying 
how faithfully it follows the Schrodinger equation. Considerable improvements in accuracy have 
been demonstrated in comparison with the usual Davydov trial states, i.e., the single Di and D 2 
Ansatze. With an increase in the number of the Davydov D 2 trial states in the multi-D 2 Ansatz , 
deviation from the exact Schrodinger dynamics is gradually diminished, leading to a numerically 
exact solution to the Schrodinger equation. 


I. INTRODUCTION 

Since the advent of the ultrafast laser spectroscopy, 
much attention has been devoted to the relaxation dy¬ 
namics of photoexcited entities, such as polarons in inor¬ 
ganic liquids and solids Eli , charge carriers in topolog¬ 
ical insulators Hi, electron and hole trapping of semi¬ 
conductor nanoparticles @-@], and electron-hole pairs in 
light-harvesting complexes of photosynthetic organisms 
|9l-[l2j. In photosynthesis, it is suggested that quantum 
coherence might play a significant role in achieving the 
remarkably high efficiency (> 95%) of the excitation en¬ 
ergy transport EMI, which can be studied by observ¬ 
ing electronic superpositions and their evolution using 
2D electronic spectroscopy techniques El HIE®- 
Emerging technological capabilities to control femtosec¬ 
ond pulse durations and down-to-one-hertz bandwidth 
resolutions provide novel probes on vibrational dynamics 
and excitation relaxation which were elusive in the past 
El- Recently, developments in ultrafast laser physics 



Figure 1: Schematic of the Holstein model. In the molecular 
ring with N — 16 sites, the circles and oscillators denote the 
excitons and phonons, respectively, and the arrows represent 
the diagonal coupling (one on one) and off-diagonal coupling 
(one on two) between the excitons and phonons. 
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and technologies allow us to study the nonequilibrium 
carrier/exciton dynamics that was previously inaccessi¬ 
ble to traditional linear optical spectroscopy. In contrast, 
modeling of polaron dynamics has not received much de¬ 
served attention. 

The Holstein molecular crystal model has been exten¬ 
sively used to study properties of polarons in molecu¬ 
lar crystals and biological systems [19]. As an exam¬ 
ple, a molecular ring with 16 sites, each coupled with 
a phonon mode, is shown in Fig. [lj Two kinds of 
exciton-phonon interactions can be included in the Hol¬ 
stein model, namely, the diagonal coupling (marked by 
an arrow attached to a single site) as a nontrivial depen¬ 
dence of the exciton site energies on the lattice coordi¬ 
nates, and the off-diagonal coupling (marked by an ar¬ 
row attached to two nearest-neighboring sites) as a non¬ 
trivial dependence of the exciton transfer integral on the 
lattice coordinates [20|. Simultaneous presence of diag¬ 
onal and off-diagonal coupling seems crucial to charac¬ 
terize solid-state excimers, where a variety of experimen¬ 
tal and theoretical considerations imply a strong depen¬ 
dence of electronic tunneling upon certain coordinated 
distortions of neighborin g m olecules in the formation of 
bound excited states [2l|, [ 23 . However, in the literature, 
little attention has been paid to the Hamiltonians con¬ 
taining the off-diagonal exciton-phonon coupling due to 
inherent difficulties to obtain reliable solutions [23|, es¬ 
pecially for the polaron dynamics [24| . Early treatments 
of off-diagonal coupling include the Munn-Silbey theory 
[25l Hfj which is based upon a perturbative approach 
with added constraints on canonical transformation co¬ 
efficients determined by a self-consistency equation. The 
global-local (GL )Ansatz EM, formulated by Zhao et 
al. in the early 1990s, was later employed in combina¬ 
tion with the dynamic coherent potential approximation 
(with the Hartree approximation) to arrive at a state- 
of-the-art ground-state wave function as well as higher 
eigenstates [29^ . 

In the absence of an exact solution, various numeri¬ 
cal approaches were developed in the past few decades, 
including the exact diagonalization (ED) [30, $1], quan- 
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turn Monte Carlo (QMC) simulation 


varia- 


onalization (YEP) (42l. kL,j. «,x*v* ^ 

coherent states [39|. Most of these approaches 


tional method [28|, l35l439| . density matrix renormaliza¬ 
tion group fPMRG ) I40l. l4ll . the variational exact diag- 

and the method of relevant 
were de¬ 
signed to probe the ground-state properties. For excited- 
state properties and dynamics of the polaronic systems, 
however, few of them provide a satisfactory resolution. 
For example, a time-dependent variant of PMRG, i.e., 
t-PMRG [434461]. was developed to elucidate the polaron 
dynamics. Yet, it cannot accurately simulate the sys¬ 
tem dynamics from an arbitrary initial state, since high- 
lying excited states can not be adequately described by 
PMRG. Fortunately, the variational approach is still ef¬ 
fective in dealing with polaron dynamics so long as a 
proper trial wave function is chosen. Previously, static 
properties of the Holstein polaron have been examined 
using a series of trial wave functions based upon phonon 
coherent states, such as the Toyozawa Ansatz [271. 1471. I48] . 


the GL Ansatz 


and a delocalized form 


of the Pavydov Pi Ansatz [50|. By using these Ansatze , 
the ground state band and the self-trapping phenomenon 
were adequately investigated. To simulate the time evo¬ 
lution of the Holstein polaron, at least two Pavydov 
Ansatze, namely, the Pi and P 2 Ansatze [5ll453j. were 
used following the Pirac-Frenkel variation scheme [54| , a 
powerful apparatus to reveal accurate dynamics of quan¬ 
tum many-body systems. Time-dependent variational 
parameters which specify the trial state are obtained 
from solving a set of coupled differential equations gen¬ 
erated by the Lagrangian formalism of the Pirac Frenkel 
variation. Validity of this approach is carefully checked 
by quantifying how faithfully the trial state follows the 
Schrodinger equation [24], [53, [56| . Numerical results show 
that the Pi Ansatz is effective and accurate in studying 
the Holstein polaron dynamics with the diagonal cou¬ 
pling, but fails to describe the case with the off-diagonal 
coupling. Pespite being a simplified version of the Pi 
trial state, however, the P 2 Ansatz instead can deal with 
the off-diagonal coupling case, albeit with non-negligible 
deviation from the exact solution to the Schrodinger dy¬ 
namics [24| . The purpose of this paper is to test the feasi¬ 
bility of using a superposition of the Pavydov trial states 
to study the dynamics of the Holstein model with simul¬ 
taneous diagonal and off-diagonal coupling. For simplic¬ 
ity, only the multi-P 2 Ansatz is probed in this work, and 
the validity of the trial state will be comprehensively in¬ 
vestigated. 


The paper is organized as follows. In Sec. II we in¬ 
troduce the model Hamiltonian and the multi-P 2 Ansatz 
for the studies of the polaron dynamics. The quantita¬ 
tive measure of the validity of the variational method is 
explained. In Sec. Ill, numerical results from our inves¬ 
tigation on the dynamics of the Holstein polaron in the 
diagonal and off-diagonal coupling regimes are displayed 
and discussed. Conclusions are drawn in Sec. IV. 


II. METHOPOLOGY 

The one-dimensional Holstein molecular crystal model 
for the exciton-phonon system can be described by the 
Hamiltonian below 


H = H ex + H p h + 4 d ‘- S P h + (1) 

where H ex , H p h and H e x-p h correspond to the exciton 
Hamiltonian, bath (phonon) Hamiltonian and exciton- 
phonon coupling Hamiltonian defined as 

Hex — J ^ ^ h^(h n _|_i T <2 n _i), 

n 

Hph = ^ ^ ^qb\b q , 

Q 

4t- S ph = + t )? 

n,q 

-ffex-ph = \<t> FA { &f n d 'n+i[e iqn (e lg - l)b q + H.C.] 

n,q 

+ata n _ 1 [e i ^(l - e~*)b q + H.c.]} , (2) 

where H.c. denotes the Hermitian conjugate, uo q is the 
phonon frequency at the momentum g, a) n ( a n ) is the 
exciton creation (annihilation) operator for the n-th 
molecule, and ( b q ) is the creation (annihilation) op¬ 
erator of a phonon with the momentum g, 

bl = N~ 1/2 J2e iqn bl St =N~V 2 Y. e ~ iq % ( 3 ) 

n q 

The parameters J,g and p represent the transfer inte¬ 
gral, diagonal and off-diagonal coupling strengths, re¬ 
spectively, and N = 16 is the number of sites in the Hol¬ 
stein polaron. In this paper, a linear dispersion phonon 
band is assumed, 


ujq = w 0 [1 + W{2\q\/n) - 1], (4) 


where cjo denotes the central energy of the phonon band, 
W is the band width between 0 and 1 , and the momen¬ 
tum is set to be q = 27 rl/N with (/ = —y + 1,..., A). 

A trial state, termed as the “Pavydov multi-P 2 
Ansatz ,” is adopted 


|D 2 M (i)> = 

N M 

ipi,n(t)exp 

n= 1 i= 1 


(5) 





ph? 


where a^(a n ) is the creation (annihilation) operator of 
a exciton at the n-th site, W q (b q ) is the creation (anni¬ 
hilation) operator of a phonon with momentum g, and 
the variational parameters ^ j7l and A ^ q denote the ex¬ 
citon probability and phonon displacement, respectively. 
Moreover, n and i represent the ranks of the site in the 
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(b) 


Figure 2: Schematics of the initial states of the Holstein po- 
laron are given in (a) and (b) for the diagonal and off-diagonal 
coupling cases, respectively. The exciton is prepared in (a) 
only at the site n = N/ 2, i.e., = 0) = N / 2 i while in 

(b) the exciton is created at the two nearest neighboring sites, 
i.e., 'i/jn = N /2 + S n ,N/ 2 + 1 )/V%- The phonon displacement 
coefficient X n (t = 0) =0 is set. 

molecular ring and the coherent superposition state, re¬ 
spectively. 

The equations of the motion are derived for the vari¬ 
ational parameters ?/An and Aby adopting the Dirac- 
Frenkel variational method, in which the Lagrangian L 
is formulated as 

L = < D 2 M C0ly^-tf|Df(i)) 

-it 

(D"(t)|^|D"(t)) - <D 2 M (i)|X|D 2 M (i)> 

- (D?(t)\H\D?(t)). 


As shown in Fig. [ 2 j the initial states of the Holstein 
polaron are prepared for the diagonal coupling case with 
the coupling strength g > </> = 0 in (a) and off-diagonal 
coupling case with the coupling strength > g = 0 in 
(b). In order to avoid the singularity, a little noise uni¬ 
formly distributed from [—£,£] is added with 5 = 10 -5 
in the initial state for both ^ j7l and A^. For each set 
of the coefficients W,g,J and defined in Eqs. m and 

more than 100 initial states are used in the simula¬ 
tions. A power-law relation between the CPU time and 
the multiplicity M has been found in our time-dependent 
variational approach, and the value of the exponent 2.8 
indicates that the time complexity of the program is 
0(n 3 ). Even for the largest value of the multiplicity 
in our paper, i.e., M = 64, the memory usage of the 
multi-D 2 Ansatz is 256 MB (million bytes), still bearable 
for the computation. Though both the CPU time and 
memory usage of the multi-D 2 Ansatz are much larger 
than those of the single D 2 Ansatz (less than 1 hour 
for the CPU time and 0.1 MB for the memory usage), 
the multi-E^ Ansatz has substantially improved the ac¬ 
curacy on variational dynamics. The energy of the sys¬ 
tem -E/totai — E/gx + E p h + -Edi a g H - E 0 fi is calculated based 
on the multi-D 2 Ansatz in Eq. m, 

E ex = <Df|i? ex |Df> 

M 

= - J EE ^f,nW'*,n+1 +1pi,n-l)Sj,i, 

i,j n 

E ph = (Df|4 h |Df) 

M 

ij n q 

M 

Ediag = (D 2 M |4t a 4hl D 2 M ) = -^EE€n^n 

hJ n 

^^(e^A^ + e-^A*J%, 
q 

M 

E oS = (D 2 M |iJ- d j ph |D 2 M ) = -^^^ W9 % i 

i,j n 

{^l n ^ n+1 [e iqn (e iq -l) A^+H.c.] 
+P jtn il>i, n -i[e iqn (l - e~ iq )Xi tq + H.c.]} , 

where Sj^ = (Aj|A^) is the Debye-Waller factor defined 
as 


ih 

~2 


From this Lagrangian, the equation of the motion for a 
and its time derivative d(t) can be obtained, 


d i dL\ dL 

dt \dd*J da * 


(7) 


where a is one of the variational parameters and 
A i :q in Eq. (0. Details on derivation of the equations 
of the motion for the Holstein polaron dynamics with a 
multi-D 2 Ansatz are given in Appendix A. 



The normalization of the wave function Norm = 
(D^ID^) is also calculated for conservation. 

Furthermore, the exciton probability P ex (t,n) and the 
phonon displacement X p h(t, n) are also calculated for the 
dynamics of the Holstein polaron. Firstly, the reduced 
single-exciton density matrix p mn (t) = Trfp^a^an] is 








obtained by solving the coupled equations of variational 
parameters, where p(t) = |D^)(D^| is the full density 
matrix at the zero temperature. After substituting the 
trial state (t)) of Eq. (j5j), the reduced density matrix 
is then derived as 


M 

Pmn(t ) = 'tpj , mnSj,i ■ (10) 

hi 


Thus, the exciton probabilities P ex (t,n) = p nn (t) (n = 
1,2,..., TV) are obtained from the diagonal elements of 
the reduced density matrix. The phonon displacement 
A^ph (£, n') in real space at the n'—th site is calculated by 


x ph (t,n') = (H) 


-L F e i9n ' 

-L W e -w»' 


«,j n 


E E 


i,j n 


Optical spectroscopy is also an important aspect of the 
polaron dynamics, as it can provide valuable information 
on various correlation functions. In this work, the linear 
absorption spectra for the polaron dynamics calculated 
with different Ansatze have been studied to check the va¬ 
lidity of these trial wave functions. The linear absorption 
spectra F(w) can be obtained by the Fourier transforma¬ 
tion, 



Figure 3: F ex , F p h, Fdi ag and E to tai obtained by the multi- 
D 2 Ansatz with M — 16, are displayed as a function of the 
time t for a molecular ring with N = 16 sites. The param¬ 
eters including the transfer integral J — 0.1, phonon energy 
bandwidth W = 0.5, diagonal coupling strength g — 1 and 
off-diagonal coupling strength (f> = 0 are set. 


where |D ^{t)) and |T(£)) obey Eq. (J7j) and the 
Schrodinger equation d\^(t))/dt = ^JT|T(t)}, respec¬ 
tively. Using the Schrodinger equation and the relation¬ 
ship |T(t)) = |Df (£)} at the moment £, the deviation 
vector | S(t)) can be calculated as 



where F(t) is the autocorrelation function of the exciton- 
phonon system, which is defined as 

F(t) = ph(0|ex(0|e i - &t Fe-^ t pt| 0 ) ex |0) ph 

= ph(0|ex(0| J Pe-^ t J pt| 0 ) ex |0) ph , (13) 


|J(t)) = ^|D"(t))-^|D"(t)). (16) 

Thus, deviation from the exact Schrodinger dynamics can 
be indicat ed by the amplitude of the deviation vector 
A (t) = (<5(£)|5(t)). In order to view the deviation in 

the parameter space (IF, J, g, 0), a dimensionless relative 
deviation a is calculated as 


with the polarization operator 

P = H l°)ex ex(0| + |0) ex e*(0|ap. (14) 

n 

Details on how to calculate the linear absorption spec¬ 
tra of a one-dimensional exciton-phonon system from the 
multi-D 2 , single D 2 and single Di Ansatze are given in 
the Appendix C. 

Finally, the validity of our Ansatz , Eq. (5), will be 
closely scrutinized. Assuming the trial wavefunction 
ID 2 d (t)) = |T(£)) at the time £, a deviation vector | S(t)) 
is then introduced to quantify the accuracy of the varia¬ 
tional method, 


max{A(t)} 
mean{£ , ph (t)} ’ 


t G [0, £ max ]. 


(17) 


where the phonon energy E p h(t ) is the main energy of 
the system, and is almost the same with the amplitude 
of the time derivative of the wave function, 


N err (t) = 


Ud^)|£ 2 | D m W) 

A E, 


( 18 ) 


l^)) = |l^W)-||DfW). 


(15) since (E) = (D^ 1 (t)\H(t)\D^(t)) « 0 in this paper. 
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III. NUMERICAL RESULTS 


A. Diagonal coupling 


The long-time behavior of the Holstein polaron dynam¬ 
ics is described by Eq. ©• Fig. [3] shows the evolution of 
the system energies, including the total energy E t otai, 
the phonon energy Eph, the exciton energy E ex , and the 
exciton-phonon interaction energy E^ ag , for the diago¬ 
nal coupling case with J = 0.1, W = 0.5, and g = 1. For 
N = 16 sites in the molecular ring, the Ansatz is formed 
from superposition of M = 16 D 2 wave functions, and 
the initial state as shown in Fig. [2fa) is used. The peri¬ 
odicity of the system energies is given by T = Aujo/tt, in 
perfect agreement with the expectation of N/4W. The 
finding that Ep^ ~ E^[ ag and E eyi — E^ 0 ^ a \ ~ 0 shows 
that the total energy is conserved in the Dirac-Frenkel 
variational dynamics. 

The exciton probability P ex (t,n) and the phonon dis¬ 
placement X p h(£,n) from the multi-D 2 Ansatz are com¬ 
pared with those from the single D 2 and Di Ansatze. The 
latter can be written as 


|Di(t)> = 

N 

P't/yft)ffln|0)exexp 
n=1 


(^n,qb\ 


(19) 


|0) p h, 


where the displacement coefficient \ n , q is not only depen¬ 
dent on the moment g, but also on the site index n in the 
molecular ring. As referred in the “Introduction,” while 
the Di Ansatz is effective in the diagonal coupling case, 
it fails to describe the polaron dynamics of the Holstein 
model with off-diagonal coupling. 

In Fig. 21 the time evolution of the exciton probabil¬ 
ity P ex (t,n) and the phonon displacement X p h(t,n) are 
shown in a weak-coupling case with g = 0.1 calculated 
using the D 2 Ansatz [panels (a) and (b)], the D^^ 16 
Ansatz [panels (c) and (d)], and the Di Ansatz [panels 
(e) and (f)]. Similarly, time-dependent behaviors of both 
P e x{t,n) and A p h(t,n) are found to be almost the same 
in the multi-D 2 and the single Di Ansatze , but at vari¬ 
ance with those in the single D 2 Ansatz. Moreover, the 
propagation of the exciton wave packets can be found in 
Fig. H c) with the velocity v ~ cd>o /27r, consistent with 
that obtained by the Di Ansatz in Fig. 0(e). 

The behavior of the Holstein polaron in the interme¬ 
diate diagonal coupling regime is shown in Fig. [5] at 
the diagonal coupling strength g = 0.2. Similar time- 
dependent behavior in P ex (t,n) and A p h(t,n) is spotted 
in the multi-D 2 Ansatz and the Di Ansatz , but not in 
the single D 2 Ansatz. The speed of the localized exci¬ 
ton wave packets v ~ ujq/Att is then measured from both 
Fig. [5] (c) and (e) to be nearly half of that in the weak 
coupling case of g = 0.1. It suggests that the veloc¬ 
ity v is inversely proportional to the diagonal coupling 
strength g. Moreover, similar phonon propagation pat¬ 
terns are found in (d) and (e), including the sound waves 




sites sites 


Figure 4: Time evolution of the exciton probability P ex (t,n) 
and the phonon displacement X p h(t,n) obtained by the D 2 , 
I 0 and Di Ansatze are displayed in (a)-(f) for a weak 
coupling case of J — 0.1, W — 0.5, g — 0.1 and 0 = 0. 


and the movement induced by the exciton-phonon inter¬ 
action. The latter is found to be with the same velocity 
as that of the exciton wave packets. 

The results of the variational dynamics in the strong 
coupling case with g = 0.4 are shown in Figs. Eta)-(f), 
corresponding to the D 2 , D^ -16 and Di Ansatze , respec¬ 
tively. Different from the weak and intermediate coupling 
cases in Figs. 0] and [5j all of the exciton probabilities 
Pex{t,n) and phonon displacements A p h(t,n) obtained 
by these three kinds of the trial wave functions nearly 
identical, indicating that all of the D 2 , multi-D 2 and Di 
Ansatze are effective in the strong diagonal coupling case. 

Using the amplitude of the deviation vector A (t) from 
the exact Schrodinger dynamics, the validity of the multi- 
D 2 Ansatz is investigated for various numbers of M de¬ 
fined in Eq. (|5j), as shown in Fig. [T] The amplitude 
of A(£) is almost constant, and the time-averaged value 
(A (£)) monotonically decreases with M. It indicates that 
the multi-D 2 trial state approaches the exact solution to 
the Schrodinger equation when M is increased. For com¬ 
parison, A(t) obtained by the Di Ansatz is also plotted. 
The time-average value of A (t) is smaller than that ob¬ 
tained from the single D 2 Ansatz (M = 1), consistent 
with our previous contention that Di Ansatz is more ac¬ 
curate than D 2 Ansatz in the diagonal coupling regime. 
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P Jt,n) by D. 




PJt,n) by 16-D 


X ph (t,n) by 16-D ? 



P ex (t,n)byD, 


x ph (t.n)byD, 



Figure 5: Time evolution of the exciton probability P ex (t,n) 
and the phonon displacement X p h(t,n) obtained by the D 2 , 
D^ =16 and Di trial states are displayed in (a)-(f) for J — 
0.1, W = 0.5, g = 0.2, and 0 = 0. 


Figure 6: For the strong diagonal coupling strength g = 0.4, 
time evolution of the exciton probability P ex (t,n) and the 
phonon displacement A p h (t,n) obtained by the D 2 , D^ =16 
and Di Ansatze are displayed in (a)-(f) at J — 0.1, W — 0.5 
an 1 J " 


Interestingly, (A(£)) calculated from the multi-D 2 Ansatz 
with M = 16 is 0.05, much smaller than 0.14 from the 
Di Ansatz and 0.40 from the single Di Ansatz , showing 
the superiority of the multi-D 2 Ansatz. 

Among the three trial states, the validity test of the 
multi-D 2 states is comprehensively performed for various 
values of the diagonal coupling strength g and transfer 
integral J in the diagonal coupling only regime. In Fig.[8j 
the relative deviation a defined in Eq. (ED is displayed as 
a function of the diagonal coupling strength g for various 
values of M for the case of J = 0.1, W = 0.5 and 0 = 0. 
As M increases, the relative error a decreases, especially 
for the weak coupling case of g = 0.1. For comparison, 
<7 calculated from the Di Ansatz is also shown, with a 
obviously larger than that of the multi-D 2 Ansatz with 
M = 32, further supporting the superiority of the multi- 
D 2 Ansatz over the Di state. Moreover, a ~ 0.01 at 
M m 32 indicates that the variational method based on 
the multi-D 2 Ansatz is possible to be numerically exact in 
the limit of M —>■ 00 , where “numerically exact” means 
the relative error a = 0 within numerical errors. 

The behavior of the relative error cr as a function of the 
transfer integral J is also investigated for several values of 
M. As shown in Fig.[9j a monotonically increases with J. 
When J is larger, the wavefunction approaches a plane 



Figure 7: The amplitude of the deviation vector A (t) from 
the Schrodinger equation is plotted as a function of the time 
t with the unit 2'k/ujq at J — 0.1, W = 0.5, g — 1.0 and 0 = 0. 
Lines represent the results obtained from the multi-D 2 Ansatz 
with various values of M, and the stars denote those from the 
Di Ansatz. 






































































wave, which is difficult to be described by superpositions 
of coherent states. With an increase in M, however, the 
relative deviation a is monotonically reduced indicating 
the improved efficiency of the D 2 Ansatz even in the case 
with a large transfer integral. Moreover, the results show 
that the multi-D 2 Ansatz is at least as accurate as the Di 
Ansatz if M = 32 superpositions are used. 


B. Off-diagonal coupling 

In this section, we probe the dynamics of the Holstein 
polaron in the off-diagonal coupling regime via the multi- 
D 2 Ansatz , in comparison with those obtained by the sin¬ 
gle D2 and the single Di Ansatze. The initial state takes 
the double excited states shown in Fig.[2fb). For simplic¬ 
ity, only the off-diagonal coupling is considered in the 
simulations. Time evolution of the exciton probability 
P e x(£, n) and phonon displacement X p h (£, n) obtained by 
the single Di, the single D2 and D ;^ 16 Ansatze are dis¬ 
played. In Fig. flQla) and (b) for the Di Ansatz , one can 
find that self-trapping occurs at the 8-th and 9-th sites, at 
variance with the delocalization expectation of the Hol¬ 
stein dynamics in the off-diagonal coupling case. It is 
consistent with the previous impression that Di Ansatz 
fails to describe the dynamics of the Holstein polaron in 
the off-diagonal coupling regime. 

In contrast, the spread of the exciton is found in the 
results obtained by the single D2 and multi-D 2 Ansatze. 
However, as shown in Figs. flQl c) and (d), the wave func¬ 
tion obtained by the single D2 Ansatz is still localized, 

diffhrpnt. from thnsp nhtainprl hv flip mnlti -D^ = ^ An.zn.t.7. 

st □ 



Figure 8: The relative deviation a defined in Eq. (H3 is dis¬ 
played as a function of the diagonal coupling strength g at 
J = 0.1, IT = 0.5 and 0 = 0. The lines with dots represent 
the results obtained from the multi-D 2 Ansatz with various 
values of M, and the dashed line denotes those from the Di 
Ansatz. The size of the molecular ring is N = 16. 



J 

Figure 9: The relative deviation a in a 16-site molecular ring 
is displayed as a function of the transfer integral J in the case 
of g — 1, IT = 0.5 and 0 = 0. The lines with dots represent 
the results obtained from the multi-D 2 Ansatz with various 
values of M, and the dashed line denotes those from the Di 
Ansatz. 


Pjtw, yt,n)by Dl 



P ex (t,n> by D 2 X ph (t,n)byD 2 




sites 


sites 


Figure 10: In the off-diagonal coupling case with 0=1 and 
J — W — g — 0, time evolution of the exciton probability 
Pex(t,n) and phonon displacement X p h(t,n) obtained by the 
lDi, D 2 and D^ =16 Ansatze are displayed in (a)-(f). 
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1/M 


Figure 11: The relative deviation a from the multi-D 2 Ansatz 
is displayed as a function of 1 /M for the the off-diagonal cou¬ 
pling case with the strength 0=1, and other parameters 
aJ = W = g = 0 are set. The size of the molecular ring is 
N = 16. In the inset, the relationship a ~ M~ v is displayed 
on a log-log scale, and the dashed line represents a power-law 
fit. 


multi-D 2 Ansatz is more effective than the single D 2 
Ansatz for depicting the Holstein polaron dynamics in the 
off-diagonal coupling case. Moreover, one can find three 
different stages of the exciton motion in figure HDTe), 
separated by the characteristic time t% « 2tt/ujo and 
£2 ~ Stt/ujq. When t < t\, the exciton is self-trapped in 
the initial state, pointing to the localized exciton state. 
After that, the exciton starts to spread over the molec¬ 
ular ring with the velocity v ~ ujq/tt. Until t > £ 2 , a 
uniform distribution of the exciton wave packets appears, 
indicating that the exciton is in a delocalized state. 

Via the relative deviation a defined in Eq. (H7|) . the 
validity of the multi-D 2 Ansatz can be further confirmed. 
In Fig. [lTJ the relative deviation a is plotted as a func¬ 
tion of 1 /M for the off-diagonal coupling only case with 
< f> = 1.0. Both the transfer integral J and the phonon 
bandwidth are set to be zero. As M increases, the rela¬ 
tive deviation a decreases and approaches zero as M goes 
to infinity. For example, the value of cr(M = 60) = 0.26 
is much smaller than a(M = 1 ) = 0.67. According to 
the fitting in the inset, the relationship a ~ M~ v is re¬ 
vealed with the exponent v = 0.29(1), further confirming 
the prediction a = 0 in the limit of M —>■ 00 . Hence, it 
can be concluded that the variational method based on 
the multi-D 2 Ansatz is possible to be numerically exact 
(<r = 0) in both of the diagonal and off-diagonal coupling 
regimes. Since any quantum state of a system of multiple 
oscillators (boson modes) can be represented by a con¬ 
tinuous superposition of coherent states (often referred 
to as the unit decomposition property of the aforemen¬ 


J=0.1, W=0.5 J=0.5, W=0.5 



J=0 1,W=0.5 J=0.5, W=0.5 



Figure 12: The relative deviation a is displayed for the sin¬ 
gle D 2 and multi-D 2 Ansdtze as a function of the diagonal 
coupling strength g and off-diagonal coupling strength 0. In 
(a) and (b), the exciton at t — 0 is created on a single site 
n = N/ 2, while in (c) and (d), a two-site occupied state is 
used as the initial state. 


tioned states), there is a good chance that the multi- 
D 2 Ansatz is numerically exact, i.e., provides exact re¬ 
sults in the M —>> 00 limit. However, it remains unclear 
how practical is the above statement, since the values of 
the multiplicity M, needed for the Ansatz to converge to 
the exact solution of the dynamical Schrodinger equation 
might be unrealistically large. The above question will 
be addressed elsewhere. 

Finally, accuracy of the multi-D 2 Ansatz is quantified 
for the parameter regime 0 < g,<fi < 1 , in comparison 
with that of the single-D 2 Ansatz , as shown in Fig. [ 12 j 
The influence of the excitonic initial state is also taken 
into account. Figs. [12] (a) and (b) correspond to the one- 
site occupied initial state shown in Fig.[2ja), and Figs. [ 12 ] 
(c) and (d) to the two-site occupied initial states shown 
in Fig. [2jb). For each initial state, two different values 
of the transfer integral J = 0.1 and 0.5 are used. Our 
results show that D^ =1 ^ Ansatz deviates little from the 
exact solutions of the time-dependent Schrodinger equa¬ 
tion in the whole parameter regime for both the two ini¬ 
tial states and the cases with small and large transfer 
integral, thereby further confirming the validity of the 
variational method. Moreover, the significant improve¬ 
ment of the validity for the multi-D 2 Ansatz from the 
single D 2 Ansatz is found especially in the weak diago¬ 
nal or off-diagonal coupling regimes, confirming the high 
accuracy of the Ansatz. 
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Figure 13: Linear absorption spectra of a 16-sites, one¬ 

dimensional ring of a coupled exciton-phonon system are dis¬ 
played in (a) for the single D 2 Ansatz , and in (b) for the 
A nsa tz. Two kinds of the initial states including the 
one-site occupied and Gaussian occupied exciton distributions 
are used in the simulations, and the Huang-Rhys factor is 
S = 2.56. For comparison, the numerical results obtained by 
the single Di Ansatz and the fitting of a Poisson distribution 
are given in (b) with the solid circles and bars, respectively. 
A rescaled factor is applied to normalize the spectral maxima 
for facilitate comparisons. 


C. Absorption spectra 

Besides the relative deviation, the validity of the 
Ansatze in providing reliable dynamical information can 
also be gauged by the accuracy of optical spectra, as an¬ 
alytical expressions for the absorption and fluorescence 
spectra are well-known if the transfer integral J and the 
phonon bandwidth W are negligible. In this study, the 
set of the parameters J = 0.1, W = 0.1, g = 0.4 and 


(j) = 0 is used, and the Huang-Rhys factor S is then cal¬ 
culated according to the relaxation energy defined by 

Er = f°° AAA du = J2 fffa = Sujo, (20) 

J-00 w q 

where ujo is the central energy of the phonon band, 
9q = 9 = 0.4 is the diagonal coupling, and uj q is the 
frequency at the moment q. From this equation, we 
can obtain S = 2.56 corresponding to the diagonal cou¬ 
pling strength g = 0.4. Two types of initial states, 
one-site-occupied excited state and the excitation with a 
Gaussian-type distribution spanned on 7 sites, have been 
applied in the study of optical spectra. To facilitate com¬ 
parisons, spectral maxima are normalized to unity. 

Three trial states including the single Di, the single 
D 2 and the multi-D 2 Ansatze are investigated, as they 
differ in terms of the variational parameters in describ¬ 
ing the phonon behavior. Linear absorption spectrum is 
a very useful indicator of the Ansatz validity in the in¬ 
vestigation of the dynamics of a polaron system. For the 
single D 2 Ansatz , the phonon displacement is only de¬ 
scribed by one set of variational parameters, leading to 
the lack of exciton-phonon correlation between exciton 
and phonon, and eventually to the absence of long-range 
correlation in the autocorrelation function and an inac¬ 
curate description of optical spectra. Only for the cases 
with the one-site occupied initial state under strong diag¬ 
onal coupling and small J where the exciton is localized, 
as shown by the black solid line in Fig. H5T a). the in¬ 
ability of the single D 2 Ansatz is alleviated. As for the 
red line, where the initial electronic excitation adopts a 
Gaussian distribution, the spectrum even exhibits nega¬ 
tive values around uj = — 3.5co>o- In contrast, the single 
Di and the multi-D 2 Ansatze guarantee the long-range 
exciton-phonon correlation, therefore can provide accu¬ 
rate absorption spectra. Fig. H5Tb) shows similar correct 
spectra for both single Di and multi-D 2 Ansatze. 

According to the Huang-Rhys theory, the phonon 
sidesbands at zero temperature follow a Poisson distri¬ 
bution, 

00 on 

F (uj) = exp (— S) — -S (uj + E r — nuj) (21) 

n n - 
n =0 

From Eq. (l2l|h the left most sideband, n = 0, corre¬ 
sponding to the zero-phonon line, should be located at 
uj = —Sujo, ie. uj = — 2.56co>o consistent with our result 
uj = — 2.64co>o as shown in Fig. H5T b). What is more, the 
tallest peaks at n = 1 and 2 show almost same height, in 
agreement with the predict that tallest of phonon side¬ 
bands should ben = S —1 = 1.56 peak when S 1 . Fur¬ 
ther, using a fitting method, the bar plot of the Poisson 
distribution with the parameter A = 2.2 is given, con¬ 
sistent with our spectra obtained from time-dependent 
evolution of single Di and multi-D 2 Ansatze. The slight 
deviation of the Poisson parameter A from the Huang- 
Rhys factor S = 2.56 is due to the nonzero values of W 
and J. 
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IV. CONCLUSION 

Using the multi-D 2 Ansatz as the trial state, we have 
systematically studied the dynamics of a one-dimensional 
Holstein polaron with simultaneous diagonal and off- 
diagonal exciton-phonon coupling via the Dirac-Frenkel 
time-dependent variational approach. Compared to the 
single D 2 Ansatz , the multi-D 2 counterpart is much more 
sophisticated and contains more flexible variational pa¬ 
rameters, leading to superior quality simulations of po¬ 
laron dynamics. Special attention is paid to testing 
the validity of our time-dependent variational approach 
by quantifying how closely the trial state follows the 
Schrodinger equation. Linear absorption spectra derived 
from the trial state are also studied as a sensitive indica¬ 
tor of the Ansatz validity in the investigation of polaron 
dynamics. 

Our numerical results show considerable improvements 
in accuracy of polaron dynamics by the multi-D 2 Ansatz , 
in comparison with the usual, single Di and D 2 trial 
states. In the diagonal coupling regime, the multi-D 2 
Ansatz is found to be at least as accurate as the sin¬ 
gle Di Ansatz for various values of the diagonal-coupling 
strength g and the transfer integral J, and remarkably 
better than the single D 2 Ansatz. In the off-diagonal cou¬ 
pling regime, however, the multi-D 2 Ansatz is shown to 
be much more potent in depicting the Holstein polaron 
dynamics than the single D 2 Ansatz , while the single Di 
Ansatz fails completely. As the number of superposed 
states M increases, one can find visible decays of the rel¬ 
ative deviation a in the weak diagonal coupling regime 
as well as the off-diagonal coupling regime, confirming 
respectable accuracies of the multi-D 2 Ansatze. More¬ 
over, a = 0 is predicted by the numerical fitting in the 
limit of M —y oo, inferring that the Dirac-Frenkel time- 
dependent variational approach based on the multi-D 2 is 
possible to be numerically exact. 

The single Davydov Di Ansatz is a trial state with suf¬ 
ficient flexibilities to handle accurately quantum dynam¬ 
ics from the celebrated spin-boson model to larg e light¬ 
harvesting complexes in photosynthesis [57], l58j . Very 
recently, a systematic coherent-state expansion of the 
ground state wave function that is based on the Davydov 
Di Ansatz , which we shall call the “multi-Di Ansatz ,” is 
developed for a number of models [591461) . It is a general¬ 
ization of a variational wave function originally proposed 
by Silbey and Harris [HI, and also an extension of the 
hierarchy of translationlly-invariant Ansatze proposed by 
Zhao et al. 00. The results of the quantum phase 
transition obtained from the multi-Di Ansatz are more 
accurate than that from the single Di Ansatz , and they 
are in agreement with DMRG and ED results, showing 
the superiority of the multi-Di Ansatz. The successful 
application of the Multi-D 2 Ansatz in the Holstein po¬ 
laron dynamics points to the possibility that the multi-Di 
Ansatz is not only valid for studying static properties of 
the model Hamiltonians, but also holds promise to their 
dynamics simulation. Our work here therefore serves as 


a proof of concept. Dynamics simulation of a Holstein 
polaron by the multi-D 2 Ansatz has convincingly shown 
that even a relatively simple wave function such as the 
Davydov D 2 Ansatz , when used in an expandable super¬ 
position, can still produce superior results. If the D 2 
Ansatz were to be replaced by the more sophisticated Di 
trial state, and applied in a multitude as demonstrated 
here, much better results can be expected on simulating 
quantum dynamics of many-body systems. Work in this 
direction is now in progress. 
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Appendix A: Time evolution of the multi-D 2 trial 
state 


As mentioned in Sec.II “Methodology”, the time evolu¬ 
tion for the multi-D 2 Ansatz can be derived by employing 
Dirac-Frenkel time-dependent variational method. Ac¬ 
cording to the definition of the multi-D 2 Ansatz in Eq. ([5]) 
and the Dirac-Frenkel variational principle, the varia¬ 
tional parameters ^ ?n (£) and A( t ) should obey 


d_ ( dL \ 

dt v a dJ 

d_ ( dL \ 

dt \ 9 K,q) 


dL 

Wn 

dL 


0 , 

0 , 


(Al) 


where L is the Lagrangian defined in Eq. ©. Inside the 
Lagrangian, the first term reads as following, 


(Df (t) | tL | D m ^ _ ( D M (t) | | d m 

M 

i,j n 


M 




i,j n 

1 • 


q L 


-(A* g A Jg + X* q Xj q ) 


S" ^iq\q) S" ^jq^iq 




Sji, (A2) 


and the second term is 

(D^(t)\H\D^(t)) = E ex + E ph + E diag + E off ( A3) 

where E ex , E p h, E d i ag and E 0 ^ denote the energies of the 
exciton, phonon, diagonal couple term and off-diagonal 
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coupling term, respectively, which can be calculated 
based on Eq. ©• 


Time partial derivatives of A i q (t) and ipi n {f) are 
then obtained by substituting Eqs. (ixa and (IVfli into 
Eq. m, which are 


and 


M 


-EE iq^ki 


M 


. M 


-EE EE «„V-lnA „S t ,, 


S( : 


2\k p \ip A ip\ p ^ip\p^ 


M 


— J ^ ^ ^ ^ 'tpkn (V^,n+1 T V^,n—l) ^iq$k,i (A5) 

i n 

— ^ ] y ] '^kn'&in | ^<7 T ^ ^ ^p^kp^ip J ^iq^ki 

in \ P / 

M 

+g EE ipl n ip in u q e lqn S k i 

i n 
M 

+g EE«»A i? E (e ijm A ip + e~ ipn \* kp ) S kli 

in p 

M 

~2 EE W A [^n + ie- i9 ”(e-^ - 1) 


+^n-ie- iqn (l -e iq )\s ki 

1 M 

-bEE ,71+1^,72 T V’/cnV^n+l) A^g 


[e ipn (e ip - l)A ip + e~ ipn {e~ ip - 1)A^] S k>i , 


M 


1 y ^ ^Pin 

. M 

2 X>*E ( 


2A* kq A ig - A iq A* q - A iq A * q ) S M 


M 

= z,n+1 T Tpi,n— 1 ) 5 fei (A4) 

i 

M 

- E E ^q^kq^iq^ki 

i Q 

M 

+9 E E ( ei9 ” A ^ + e “^ A D 

* <2 

1 M 

“^EE^Vv+i [e* 9 "(e* 9 - l)Ai, 

i <7 

+e -^ (e - ig _ 1)A * g ] Sfc . 

1 M 

- 2 ^EE w ^i,n-i [e* 9 "(l - e- i9 )A i9 
i q 

+e~ ign (l - e iq )A* kq ] S ki , 


where Sk,i is the Debye-Waller factor defined in Eq. ©• 
Unfortunately, both Eqs. (1A4[) and (IA5|) contain the cou¬ 
pling time partial derivatives of A i q (t) and ipi n (t) with 
the form JA 4 ^™ + M\ q + J 2 i, q A ^K q = B, where 
A 1 , A 2 , A 3 and B are coefficient vectors, irrelative to the 
time partial derivatives of the variational parameters. By 
numerically solving these linear equations at each time £, 
one can calculate the values of ipi n and A i q accurately. 
Runge-Kutta 4-th order method is then adopted for the 
time evolution of the Holstein polaron, including the en¬ 
ergies E/totai(^) — Eex ~b "^ph "b -Ediag T T/ 0 ff, the normal¬ 
ization Norm{t ), the exciton probability P ex (t,n) and 
phonon displacement X p h(t,n). 

Finally, the amplitude of the deviation vector A (t) de¬ 
fined in Eq. (USD is calculated as 


A 2 (t) = (S(t)\S(t)) 


M 


= EE 


n i,j L q 


E( A* nq + B* jnq A iq ) 


Jjl 


(A 6 ) 


y y {Ainq T Bi nq \j q ) 


M 


BjnqSjiBinq , 


n,g 


where the matrixes Aj nq and Bj nq are respectively de¬ 
fined as 







Figure 14: In the diagonal coupling case, the total energy of 
the system E to tai and the normalization of the wave function 
Norm are plotted as a function of the time t for a molec¬ 
ular ring with N — 16 sites. The time unit 2tt/uo denotes 
the vibrational period of the phonon. The parameters includ¬ 
ing the transfer integral J = 0.1, phonon energy bandwidth 
W — 0.5, diagonal coupling strength g — 1 and off-diagonal 
coupling strength 0 = 0 are set. 



t / (27l/CO 0 ) 


Figure 15: In the off-diagonal coupling case with 0=1 
and J — W — g — 0, the phonon energy E p h, off-diagonal 
exciton-phonon interaction energy E 0 r and total energy E tota \ 
obtained by the multi-D 2 Ansatz with M — 16, are displayed 
as a function of the time t for a molecular ring with N = 16 
sites. In the inset, the normalization of the wave function 
Norm is plotted. 


and 


A. 


jnq 

(t) 

1 


i cy^j.n (t) \j,q (t) ^j,q W + ^j,q W ^j,q (t) 


N i' l l ; 3 : n +1 W (0] 


+ 91 pj,n (t) UJ q e lQn Xj, q 

-\<M> 3 ,n +1 (t)u q e iqn (e iq - 1) A jtQ (t) 

( t ) u q e lqn (l - e~ lq ) \ jt g ( t ), 

Bjnq 

j,n (t) ^ j,q (t) 

^Pj,n {t) ^q^j,q 
+g^j,n (t) 0 J q e~ iqn 

- l -Ui,n+ 1 (t)co q e- iqn (e~ iq — l) 


the wave function Norm(t) are plotted as a function of 
the time t in Fig.[l4]for the precision test of the multi-D 2 
Ansatz with M = 16. One can find that the deviations 
of E t otai and Norm from the initial values are negligi¬ 
bly small, suggesting that the numerical results obtained 
by the Dirac-Frenkel variational dynamics based on the 
multi-D 2 trial states are reliable. 

Besides, the dynamic behavior of the Holstein polaron 
for the off-diagional coupling case is also investigated at 
J = 0,VF = 0,g = 0 and 0 = 1. As shown in Fig. [15l 
aperiodic behaviors of the system energies are found, 
consistent with the prediction of the long period time 
T —ex) due to the vanishing of the band width W = 0. 
^totai(t) = £0h + £ 0 ff ~ 0 shows the system total energy 
is conserved. In the inset, the normalization Norm(t) is 
also displayed for the conservativeness test. 


Appendix C: LINEAR ABSORPTION 


Combing the Eqs. (fT3l) and (Q3j), the autocorrelation 
function is derived 


Appendix B: Energy Conservation 

In the diagonal coupling case of J = 0.1, W = 0.5 and 
g = 1, the total energy E tota i(t) and the normalization of 


fW=f 2 EE ph <°lex <0| |0) ex |0) ph 

n m 

(Cl) 
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Using the periodic condition of the Hamiltonian H , one 
has 

ph (0|ex (0| |0) ex l^)ph 

m 

= (0| ex (0| ^ m-n e % 1 ^n-n |0) ex l^)ph 

m 

= E Ph (°lex <0| & m e~ iftt al |0) ex |0) ph . (C2) 

m 

Substituting Eq. (1C2[) into Eq. (|C1|) . one can obtain 
F (t) = E P h (0| ex <0| a m e- i6t al |0) ex |0) ph , (C3) 

m 

where e _2i ^do |0) ex |0) h is the time evolution of wave 

function from the initial state dj |0} ex |0) ph , which can 
be approximated by a Davydov trial state, for example, 
by the Multi-D2 trial state, 

e~ im a\ |0) ex |0) ph 

M 


For the single trial state, the time evolution of wave 
function from the initial state can be approximated to 

0 -iHt~ f 


e~ lHt a' 0 |0) ex |0) ph 

~ E (*) & n 


exp 


E [ x i (*) b l 


-H.c. 


|0) e JO) ph , (C5) 


which leads to the autocorrelation 


EE^.» 


exp \ Y, [\ q (t) SJ - H.c] 1 |0) ex |0) ph . (C4) e * % |0) ex l°> P h 


Ph <oi 

n m 

exp |E [ A « W h l - H c ■ | l°)ph 
= i-> 2 N E r 'i ,m W ex pE i a 9 wi 2 } 

m l q J 


Finally, the single D\ trial state is used to calculate 
the time evolution of the wave function 

-iHt* t 


The autocorrelation F (t) of the multi-Z)2 Ansatz is then 
calculated by substituting Eq. (1C4[) into Eq. (|C3 [h 


M 


^w = k 2 EEE^w P h <°i 

ij n m 

exp{E[y. {t)b q -H.c] | 

exp |E [ Ai -« (*) ^9 _ Hx ‘ | l°)ph 

M 

= e E w 

ij n 

E [ _ (I A J<7 1 2 + l A «? I 2 ) + A iq A ' 


'y ] (t) 1 

n 

exp |e; [An, (t) St - ff.c.] | |0) ex |0) ph , (C6) 

and the autocorrelation F(t) is then obtained by 
F (t) = 'y ^ y ^ {t) ph (0| 

n m 

exp |E[ A -9W^-^ C - |l°)ph 
= M 2 ^E^m Wexp ( -^El Am « Wl 2 ) • 

m \ q J 


exp 
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